\documentclass[10pt, a4paper, titlepage, oneside]{article}
\usepackage[top=2cm, bottom=1.5cm, left=2cm, right=1.5cm]{geometry}
\usepackage{graphicx}
\usepackage{pgffor}
\usepackage{amsmath}
\usepackage{mathtools}
\usepackage{textcomp}
\usepackage{multirow}
\usepackage{lscape}
\usepackage{listings}

% dot
%\usepackage{dot2texi}
%\usepackage{tikz}
%\usetikzlibrary{shapes,arrows}

% Font, accents
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{times}                % Otros: palatino, bookman, times

% Interactive index
\usepackage{hyperref}

% Spanish titles
\usepackage[spanish]{babel}

\lstset{language=Matlab}

% Make lists without bullets
%\renewenvironment{itemize}{
% \begin{list}{}{
%  \setlength{\leftmargin}{1.5em}
% }
%}{
% \end{list}
%}

% Document properties
\title{75.12 - Análisis Numérico - Grupo 2 - TP2}
\author{Grupo 2\\
	Federico Churca-Torrusio\\
	Emanuel Cruz\\
	Miguel Serigos}
\date{2013 - Verano\\
	Fecha de entrega: 03 de marzo de 2013}
\hypersetup{
  colorlinks = true,
  urlcolor = red,
  pdflang = es,
  pdfauthor = {Grupo 2 - Federico Churca-Torrusio, Emanuel Cruz, Miguel Serigos},
  pdfproducer = {Federico Churca-Torrusio, Emanuel Cruz, Miguel Serigos},
  pdfcreator = pdflatex,
  pdftitle = {75.12 - Análisis Numérico - Grupo 2 - TP2},
  pdfsubject = {},
  pdfpagemode = UseNone
}

\providecommand{\abs}[1]{\left\lvert#1\right\rvert}
\providecommand{\norm}[1]{\left\lVert#1\right\rVert}

\let\stdsection\section
\renewcommand\section{\newpage\stdsection}

\begin{document}
\maketitle

%\thispagestyle{empty}
%\cleardoublepage
\pagenumbering{roman}
\setcounter{page}{1}
\tableofcontents
%\newpage
\listoffigures
\listoftables
\newpage

\pagenumbering{arabic}
\setcounter{page}{1}

\section{Introducción}

Analizaremos la variación del stock de un almacén resolviendo la EDO partiendo del siguiente modelo:
\begin{align}
	\dfrac{dx}{dt} = T_c \cdot{x} - T_v \cdot{x}
\end{align}
Donde $T_c$ es la tasa de compras y $T_v$ la tasa de ventas.

Vamos a resolver la EDO en forma analítica y en forma numérica.
Para la forma numérica usaremos resolvedores propios para los métodos de Euler, Runge-Kutta de orden 2, y Runge-Kutta de orden 4, y compararemos los resultados con resolvedores estándar de la librería OdePkg para los métodos Runge-Kutta 2,3 y Runge-Kutta 4,5.
Compararemos todos los resultados con una implementación de la solución analítica del problema para obtener los errores.

Luego de describir la solución teórica del problema, describiremos cada método implementado, el resultor que es empleado en cada método, y su estabilidad.

Las implementaciones que realizamos no tienen fijos los parámetros del problema, sino que, al igual que la familia de algoritmos de la librería estándar (sea Matlab nativo u OdePkg para Octave), toman una función derivada arbitraria en función de $(t, x)$.
Dicha función derivada es generada en función de $(T_c, T_v)$ por otra función.
Esto nos permite hacer la segunda parte del informe sin reescribir los métodos resolvedores.

\subsection{Resolución analítica}

Partiendo del modelo teórico, luego de sacar factor común nos queda:
\begin{equation} \label{inicial}
	f_{(t, x)} =  \dfrac{dx}{dt} = (T_c-T_v) \cdot{x}	
\end{equation}

Luego resolvemos:
\begin{align}
	{dx}&=(T_c-T_v)\cdot{x}\cdot{dt}\\
	\displaystyle\int\dfrac{1}{x}dx&=\displaystyle\int(T_c-T_v)dt\\
	\ln{(x)} &= (T_c - T_v) \cdot{t} + c\\ 
	x_{(t)} &= e^{(T_c-T_v)\cdot{t}}\cdot{e^{c}}\\
	x_{(t)} &= e^{(T_c-T_v)\cdot{t}}\cdot{k}\\
	x_{(0)} &= e^{0}\cdot{k}\\
	k &= x_{(0)}
\end{align}

Finalmente: 
\begin{align}
	x_{(t)}&={e^{(T_c - T_v)\cdot{t}}} \cdot x_{(0)}
\end{align}

Reemplazando las condiciones $T_c = 0.1$, $ T_v = 0.2$ y $x_{(0)} = 100$:
\begin{align}
	x_{(t)} &= e^{-0.1 \cdot t} \cdot 100
\end{align}

\subsection{Método de Euler}

Es un método de paso entero, explícito y basado en una aproximación en adelanto para $x$:
\begin{align} \label{euler}
	x_{n+1}&= x_n + h \cdot{f_n}\\
\end{align}
Siendo $f_n= f_{(x_n,t_n)}$ la función derivada obtenida de la ecuación diferencial.
Luego reemplazamos según el modelo a resolver:
\begin{align}
	x_{n+1}&= x_n + h \cdot{({f={((T_c=0.1)-(T_v=0.2))}}\cdot x_n)} 
\end{align}
Con $h=0.1$ resulta:
\begin{align}
	x_{n+1}&= x_n \cdot{(1-0.1^2)} \\
\end{align}
Y finalmente, mediante el método de Euler $x_{n+1}$ es:
\begin{align}
	x_{n+1}&= 0.99 \cdot{x_n} 
\end{align}	

\subsection{Método de Runge-Kutta de orden 2}

Se usa para mejorar la precisión del método de Euler.
En este método se usan $q_1$ y $q_2$ como valores aproximados de la derivada de $x$ respecto de $t$.
Por lo tanto debemos calcular los incrementos $q_1$ y $q_2$, para luego usar su promedio para calcular $x_{n+1}$.
\begin{align}
	q_1&= h \cdot{f_{(t_n,x_n)}} \label{q1rk2} \\
	q_2&= h \cdot{f_{(t_n+h,x_n+q_1)}} \label{q2rk2} \\
	x_{n+1}&= x_n + \dfrac{q_1+q_2}{2} \label{xnm1rk2}
\end{align}
Reemplazando los valores en $q_1$ según la ecuación \eqref{inicial} :
\begin{align}
	q_1&= 0.1 \cdot{(0.1-0.2) \cdot{x_n}} \\
	q_1&= -0.01 \cdot{x_n} 
\end{align}
Nuevamente para $q_2$:
\begin{align}
	q_2&= 0.1 \cdot{(-0.1)} \cdot{(x_n-0.01\cdot{x_n})} \\
	q_2&= (-0.0099) \cdot{x_n}
\end{align}
Luego reemplazamos los valores de $q_1$ y $q_2$ en la expresión \eqref{xnm1rk2}
\begin{align}
	x_{n+1}&= x_n + \dfrac{(-0.01)\cdot{x_n}-0.0099 \cdot{x_n} }{2} \\
	x_{n+1}&= x_n \cdot{\dfrac{1.9801}{2}} 
\end{align}

\subsection{Método de Runge-Kutta de orden 4}

En este caso, el método de Runge-Kutta de orden 4 precisa dos valores aproximados de las derivadas más que el de orden 2, $q_3$ y $q_4$, siendo:
\begin{align}
	q_1&= h \cdot{f_{\left(t_n,x_n\right)}} \\
	q_2&= h \cdot{f_{\left(t_n+\dfrac{h}{2},x_n+\dfrac{q_1}{2}\right)}} \label{3}\\
	q_3&= h \cdot{f_{\left(t_n+\dfrac{h}{2},x_n+\dfrac{q_2}{2}\right)}} \label{eq4} \\
	q_4&= h \cdot{f_{\left(t_n+h,x_n+q_3\right)}} \\
	x_{n+1}&= x_n + \dfrac{(q_1+2\cdot{q_2}+2\cdot{q_3}+q_4)}{6}	 \label{eq5}
\end{align}
Para $q_1$, usamos \eqref{inicial}
\begin{align}
	q_1&= h \cdot{(T_c-T_v)} \cdot{x_n}
\end{align}
Reemplazando los valores de $T_c$, $T_v$ y $h$, obtenemos:
\begin{align}
	q_1 &= 0.1 \cdot{(-0.1)} \cdot{x_n} \\
	q_1 &= -0.01 \cdot{x_n}
\end{align}
Igualmente para $q_2$, de \eqref{3} :
\begin{align}
	q_2&= 0.1 \cdot{(-0.1)} \cdot{\left(x_n-\dfrac{0.01\cdot{x_n}}{2}\right)} \\
	&= x_n \cdot{0.1} \cdot{(-0.01)} \cdot{\left(1-\dfrac{0.01}{2}\right)}	\\
	q_2&= -0.00995 \cdot{x_n}
\end{align}
Calculamos ahora $q_3$ mediante la ecuación \eqref{eq4}, y reemplazando los valores de $T_c$, $T_v$ y $h$:
\begin{align}
	q_3&= 0.1 \cdot{(-0.1)} \cdot{\left(x_n+\dfrac{(-0.00995\cdot{x_n})}{2}\right)} \\
	&= x_n \cdot{(-0.01)} \cdot{\left(1-\dfrac{0.00995}{2}\right)} \\
	q_3&= -0.00995025 \cdot{x_n} 
\end{align}
Finalmente calculamos $q_4$:
\begin{align}
	q_4&= 0.1 \cdot{(-0.1)} \cdot{(x_n+(-0.00995025)\cdot{x_n})} \\
	&= x_n \cdot{(-0.01)} \cdot{(1-0.00995025)} \\
	q_4&= -0.0099004975 \cdot{x_n} 
\end{align}
Una vez calculados las aproximaciones de las derivadas, reemplazamos en \eqref{eq5} para encontrar la expresión de $x_{n+1}$:
\begin{align}
	x_{n+1}&= 0.9900498338 \cdot{x_n} \\
	x_{n+1}& = \dfrac{5.940299003}{6} \cdot{x_n} 
\end{align}

\section{Resultados}

\subsection{Gráficos}
Como en todos los casos el márgen de error es al menos en dos órdenes de magnitud menor al rango de valores, las curvas de valores se tornan indistinguibles unas de otras.
Por eso, decidimos incluir dos gráficos; uno con los valores esperados, y uno con los errores relativos.

\begin{figure}[!ht]
	\caption{Valores esperados}
	\begin{center}
		\includegraphics[scale=0.5]{resEsperados}
	\end{center}
\end{figure}

\begin{figure}[!ht]
	\caption{Errores relativos}
	\begin{center}
		\includegraphics[scale=0.5]{errores}
	\end{center}
\end{figure}

\subsection{Resultados de las resoluciones numéricas}
A continuación se adjuntan los resultados obtenidos por soluciones numéricas
$\tilde{x}_{(t)}$, junto a los valores esperados $x_{(t)}$, errores absolutos $\delta$, y
errores relativos $\epsilon$.
Los archivos fuente que se usaron para generar los resultados están en el
comprimido en el que se entregó el presente informe.

\begin{landscape}
	\input{tablas.part.tex}
\end{landscape}

\section{Estabilidad}

\subsection{Estabilidad del sistema}
El concepto de estabilidad hace referencia a la alteración de un sistema al variar sus condiciones iniciales.
Se considera que un sistema es más estable en cuanto una perturbación en sus condiciones iniciales se ve reflejada en su evolución de manera poco perceptible.

De esta forma, dado un sistema de ecuaciones diferenciales:
\begin{align}
	y_{(0)} &= y_0 \\
	y_{(x)} \prime &= f_{(x,y)}
\end{align}

Este tiene una solución única $y_{x}$.
Si tenemos otra solución a partir de variar sus condiciones iniciales, llamémosle ${y^\ast}_{(x)}$.
La llamaremos estable si cumple con las siguientes condiciones:

\begin{align}
	\exists {y^\ast}_{(x)} &: x \in [x_0,\infty) \\
	\forall x \geq x_0 &: \abs{y_{(x)} - {y^\ast}_{(x)}} \leq \varepsilon , \varepsilon \in (0, +\infty) \\
\intertext{Para el caso especial:}
	y_{(0)} &= y_0 \\
	y_{(x)}\prime &= \alpha_{(x)} \cdot {y} + f_{(x)} \\
\end{align}

La diferencia entre 2 soluciones cualesquiera está dada por la siguiente relación:

\begin{align}
	\abs{y_{(t)} - {y^\ast}_{(x)}} &= \abs{y_{0} - {y^\ast}_{0}} \cdot e ^ {\int_{t_0}^t \alpha _{(s)} ds}\\
\end{align}

Podemos decir entonces que este conjunto especial de funciones es estable si y solo si:

\begin{align}
	e ^ {\int_{t_0}^t \alpha _{(s)} ds} & \quad \text {está acotada} \\
	\lim_{t \to +\infty} e ^ {\int_{t_0}^t \alpha _{(s)} ds} &= 0	
\end{align}

En un caso particular donde $\alpha$ es constante, es condición necesaria que $\alpha \leq 0$ para que sea estable.
A partir la ecuación previa, podemos decir estonces lo siguiente:

\begin{align}
	\alpha &> 0 \quad \Rightarrow \text{Es inestable}\\
	\alpha &= 0 \quad \Rightarrow \text{Es estable}\\
	\alpha &< 0 \quad \Rightarrow \text{Es estable asintóticamente}
\end{align}

Mirando el problema planteado, podemos observar que tomando $T_c y T_v$ fijos caemos en el caso de $\alpha$ constante. Por lo tanto si queremos que el sistema sea estable, deberá cumplir con la condición $\alpha \leq 0$.

\begin{align}
	y_{(x)}\prime &= \alpha_{(x)} \cdot {y} + f_{(x)} \\
	y\prime &= (T_c-T_v)\cdot {y} \\
	\alpha_{(x)} &= T_c-T_v\\
	f_{(x)} &= 0
\end{align}

Dijimos que el sistema será estable si cumple con la condicion$\alpha \leq 0$,por tanto para nuestro caso particular, deberá cumplir con $T_c - T_v \leq 0$, y si queremos que sea asintóticamente estable, $T_v > T_c$.


\subsection{Estabilidad en los métodos numéricos}
A continuación se analizará la estabilidad de cada método numérico para resolver el sistema.
Así como el sistema en sí es más o menos estable si una variación en las condiciones iniciales se ve más o menos reflejado en la evolución del sistema, un método es más o menos estable si los errores de truncamiento cometidos alejan más o menos la solución numérica de la solución analítica.

\subsubsection{Euler}
Analizaremos la estabilidad del metodo, perturbando la ecuación de Euler \eqref{euler}.
\begin{align}
	x_{n+1} + \varepsilon_{n+1} &= x_n + \varepsilon_n + h \cdot \left( f_{(t_n,x_n + \varepsilon_n)}  \right)
\end{align}

Usando Taylor para aproximar la función $f_{t_n, x_n + \varepsilon_n}$ nos queda:
\begin{align} \label{tEuler}
	f_{(t_n, x_n + \varepsilon_n)} &= f_{(t_n,x_n)} + \frac{\partial f_{(t_n,x_n)}}{\partial x_n} \cdot \varepsilon_n
\end{align}

Reemplazando\eqref{tEuler} en \eqref{euler} nos queda:

\begin{align}
	x_{n+1} + \varepsilon_{n+1} &= x_n + \varepsilon_n + h \cdot \left( f_{(t_n,x_n)} + \frac{\partial f_{(t_n,x_n)}}{\partial x_n} \cdot \varepsilon_n \right)
\end{align}

De donde obtenemos que:

\begin{align}
	\varepsilon_{n+1} &= \varepsilon_n + h \cdot \frac{\partial f_{(t_n,x_n)}}{\partial x_n} \cdot \varepsilon_n \\
	\varepsilon_{n+1} &= \left(1 + h \cdot \frac{\partial f_{(t_n,x_n)}}{\partial x_n} \cdot \right)\cdot \varepsilon_n	
\end{align}

Obtenenmos así el factor de amplificación para el metodo de Euler:

\begin{align}
	FA_{Euler} &= \abs{1 + h \cdot \frac{\partial f_{(t_n,x_n)}}{\partial x_n}}
\end{align}

En nuestro caso, derivando $f_{(t_n,x_n)}$ respecto de $x_n$, obtenemos lo siguiente:

\begin{align}
		\frac{\partial f_{(t_n,x_n)}}{\partial x_n} &= T_c-T_v  \\
		FA_{Euler} &= \abs{1 + h \cdot (T_c-T_v)}
\end{align}

Para que el metodo sea estable, se debe cumplir que $FA \leq 1$.
Dado que $T_c - T_V \geq 0$. entonces $1 + h \cdot (T_c-T_v)$ es positivo y vale:

\begin{align}
	1 + h \cdot (T_c-T_v) &\leq 1 \\
\end{align}

De donde podemos asegurar que:

\begin{align}
	T_v &\geq T_c
\end{align}

\subsubsection{Runge-Kutta de orden 2}
Para analizar la estabilidad del método, calcularemos la expresión de $\varepsilon_{n+1}$ en función de $T_c$ y $T_v$.
Análogamente al desarrollo realizado para hallar la estabilidad para el método de Euler, deberemos comenzar perturbando la ecuación \eqref{xnm1rk2}, reemplazando las ecuaciones \eqref{q1rk2} y \eqref{q2rk2}:
\begin{align}
	q_1&= h \cdot{(T_c - T_v)} \cdot{x_n} \\
	q_2& = h \cdot{(T_c - T_v)} \cdot{\left(x_n + (T_c - T_v) \cdot{x_n} \right)}\\
	x_{n+1} &= x_n + \dfrac{h}{2} \cdot{(T_c - t_v)} \cdot{x_n} \cdot{\left(2+(T_c - T_v) \right)} \label{xestrk2}	
\end{align}

Llamamos $\varphi_{(t_n,x_n)}$ a la función:
\begin{align}
	\varphi &= \dfrac{1}{2} \cdot{(T_c - T_v)} \cdot{x_n} \cdot{\left(2+(T_c - T_v) \right)}
\end{align}
Obteniendo una ecuación similar a la resuelta para Euler:
\begin{align}
	x_{n+1} &= x_n + h \cdot{\varphi_{(t_n,x_n)}}
\end{align}
Luego, perturbando la ecuación, y aproximando por Taylor la función $\varphi_{(t_n,x_n)}$ al igual que en el método anterior, concluímos que:
\begin{align}
	\varepsilon_{n+1} = \underbrace{\left(1 + h \cdot{\dfrac{\partial{\varphi_{(t_n,x_n)}}}{\partial{x_n}}} \right)}_{F_A} \cdot{\varepsilon_n}
\end{align}

Para hallar $F_A$ calcularemos $\dfrac{\partial{\varphi_{(t_n,x_n)}}}{\partial{x_n}}$:
\begin{align}
	\dfrac{\partial{\varphi_{(t_n,x_n)}}}{\partial{x_n}} &= \dfrac{1}{2} \cdot{(T_c - T_v)} \cdot{\left(2+(T_c - T_v) \right)}
\end{align}

Reemplazamos:
\begin{align}
	F_A &= 1 + h \cdot{\dfrac{1}{2} \cdot{(T_c - T_v)} \cdot{\left( 2 + (T_c - T_v)\right)}}
\end{align}

Exigiremos $F_A \leq 1$ para que el método sea fuertemente estable:
\begin{align}
	1 + \dfrac{h}{2} \cdot{(T_c - T_v)} \cdot{\left(2+(T_c - T_v) \right)} \leq{1} \\
	(T_c - T_v) \cdot{\left(2 + (T_c - Tv) \right)} \leq{0} \\
	2 + (T_c - T_v) \leq{0} \\
	T_c \leq T_v - 2 \label{estrk2} 
\end{align}
Siendo \eqref{estrk2} la condición necesaria para que el método de Runge-Kutta de orden 2 sea fuertemente estable.

\subsubsection{Runge-Kutta de orden 4}
Para el análisis de la estabilidad del método en estudio procederemos de igual manera que en los casos anteriores.
Partiremos de la expresión de $x_{n+1}$ reemplazando los valores de $q_i$ con $i \subset{[1,4]}$ en función de $T_c$ y $T_v$.
Para simplificar las cuentas a contincuación utilizaremos las siguientes expresiones:
\begin{align}
	\lambda &= h \cdot{(T_c - T_v)} \label{lambdark4} \\
	p &= (T_c - T_v)
\end{align}

Luego, calcularemos $x_{n+1}$ en función de $q_1$:
\begin{align}
	q_1 &= h \cdot{(T_c - T_v)} \cdot{x_n} \\
	q_1 &= \lambda \cdot{x_n} \label{q1rk4}
\end{align}
$q_2$ en función de $q_1$:
\begin{align}
	q_2 &= h \cdot{(T_c - T_v)} \ cdot{\left(x_n + \dfrac{q_1}{2}\right)}\\
	&= q_1 + q_1 \cdot{\dfrac{h}{2}} \cdot{(T_c - T_v)}\\
	q_2 &= q_1 \cdot{\left(1 + \dfrac{\lambda}{2}\right)} \label{q2rk4}
\end{align}
$q_3$ en función de $q_1$:
\begin{align}
	q_3 &= h \cdot{(T_c - T_v)} \ cdot{\left(x_n + \dfrac{q_2}{2}\right)} \\
	&= q_1 + \lambda \cdot{\dfrac{q_2}{2}} \\
	&= q_1 + q_1 \cdot{(1 + \dfrac{\lambda}{2})} \cdot{\dfrac{\lambda}{2}} \\
	&= q_1 \cdot{\left(1 + \dfrac{\lambda}{2} \cdot{(1 + \dfrac{\lambda}{2})}\right)} \\
	q_3 &= q_1 \cdot{\left(1 + \dfrac{\lambda}{2} + \dfrac{\lambda^2}{4} \right)} \label{q3rk4}
\end{align}
Luego $q_4$ en función de $q_1$ es:
\begin{align}
	q_4 &= h \cdot{(T_c - T_v)} \cdot{(x_n + q_3)} \\
	&= q_1 + \lambda \cdot{q_3} \\
	&= q_1 + \lambda \cdot{q_1} \cdot{\left(1 + \dfrac{\lambda}{2} + \dfrac{\lambda^2}{2} \right)} \\
	q_4 &= q_1 \cdot{\left(1 + \lambda + \dfrac{\lambda^2}{2} + \dfrac{\lambda^3}{4}\right)} \label{q4rk4}
\end{align}
Ahora sí, reemplazando las ecuaciones \eqref{q1rk4}, \eqref{q2rk4}, \eqref{q3rk4} y \eqref{q4rk4} en \eqref{eq5}, obtenemos:
\begin{align}
	x_{n+1} &= x_n + h \cdot{p} \cdot{\dfrac{x_n}{6}} \cdot{\left(6 + 3 \cdot {\lambda} + \lambda^2 + \dfrac{\lambda^3}{4}\right)} 
\end{align}

Notamos ahora que si generamos la función $\phi_{(t_n,x_n)}$ definida por:
\begin{align}
	\phi_{(t_n,x_n)} &= p \cdot{\dfrac{x_n}{6}} \cdot{\left(6 + 3 \dot{\lambda} + \lambda^2 + \dfrac{\lambda^3}{4}\right)} 
\end{align}
Nuevamente tenemos la ecuación de la forma:
\begin{align}
	x_{n+1} = x_n + h \cdot{\phi_{(t_n,x_n)}} 
\end{align}

Por lo tanto, perturbando la ecuación y aproximando $\phi_{(t_n,x_n)}$, obtendremos la aproximación del error en cada paso:
\begin{align}
	\varepsilon_{n+1} = \left(1 + h \cdot{\dfrac{\partial{\phi_{(t_n,x_n)}}}{\partial{x_n}}} \right) \cdot{\varepsilon_n}
\end{align}

Debemos calcular $\dfrac{\partial{\phi_{(t_n,x_n)}}}{\partial{x_n}}$ para obtener el factor de amplificación para el método de Runge-Kutta de orden 4($FA_rk4$):
\begin{align}
	\dfrac{\partial{\phi_{(t_n,x_n)}}}{\partial{x_n}} &= p \cdot{\left(\dfrac{\lambda^3}{24} + \dfrac{\lambda^2}{6} + \dfrac{\lambda}{2} + 1\right)}
\end{align}

Por lo tanto:
\begin{align}
	FA_{rk4} &= 1 + h \cdot{p} \cdot{\left(\dfrac{\lambda^3}{24} + \dfrac{\lambda^2}{6} + \dfrac{\lambda}{2} + 1\right)}
\end{align}

Condicionamos $FA_{rk4}$ tal que su módulo sea menor o igual a 1 para que el método sea fuertemente estable:
\begin{align}
	1 + h \cdot{p} \cdot{\left(\dfrac{\lambda^3}{24} + \dfrac{\lambda^2}{6} + \dfrac{\lambda}{2} + 1\right)} &\leq{1} \\
	h \cdot{p} \cdot{\left(\dfrac{\lambda^3}{24} + \dfrac{\lambda^2}{6} + \dfrac{\lambda}{2} + 1\right)} &\leq{0} \\
	\dfrac{\lambda^3}{24} + \dfrac{\lambda^2}{6} + \dfrac{\lambda}{2} + 1 &\leq{0} \\
	\dfrac{\lambda^3}{12} + \dfrac{\lambda^2}{3} + \lambda + 2 &\leq{0}
\end{align}

Dicha condición se cumple para todo $\lambda \leq{-2,785294}$.
Recordando \eqref{lambdark4} obtenemos la relación:
\begin{align}
	\lambda \leq{-2,785294} \\
	h \cdot{(T_c - T_v)} \leq{-2,785294} \\
	T_c - T_v \leq{-2,785294} \\
	T_c \leq{-2,785294 + T_v} 
\end{align}
Por lo que cumpliendo dicha condición el método resulta fuertemente estable.

\section{Resultados con condiciones estables}

Para los tres métodos, los valores de $(Tc, Tv)$ que daban resultores más estables son los mismos $(Tc = 0.05, Tv = 0.9)$.
En los tres casos, el factor de amplificación es menor a $1$, lo que los hace fuertemente estables.

\foreach \s in {euler,rk2,rk4} {
	\begin{figure}[!ht]
		\caption{Factor de amplificación en función de Tc, Tv: \s}
		\begin{center}
			\includegraphics[scale=0.5]{estabilidad_\s}
		\end{center}
	\end{figure}
}

\begin{figure}[!ht]
	\caption{Valores esperados con condiciones estables}
	\begin{center}
		\includegraphics[scale=0.5]{estables_resEsperados}
	\end{center}
\end{figure}

\begin{figure}[!ht]
	\caption{Errores relativos con condiciones estables}
	\begin{center}
		\includegraphics[scale=0.5]{estables_errores}
	\end{center}
\end{figure}

\begin{landscape}
	\input{estables_tablas.part.tex}
\end{landscape}

\section{Tasa de venta independiente}

En el caso de que la tasa de venta sea constante y no dependiera del stock, la función derivada toma otra forma:

\begin{equation}
	f(x, t) = T_c \cdot x - T_v
\end{equation}

Esto lleva a otro modelo, con la siguiente solución analítica:
\begin{equation}
	x_{(t)} = \dfrac{e^{(T_c \cdot t)} \cdot (x_{(0)} \cdot T_c - T_v) + T_v}{T_c}
\end{equation}

Este sistema teóricamente es inestable si $T_c \neq 0$ (no tiene sentido una tasa de compra negativa):

\begin{align}
	x \prime &= a \cdot x + g_{(t)}\\
	g_{(t)} &= T_v\\
	a &= T_c \geq 0
\end{align}

Para que el sistema sea estable, tiene que darse $T_c = 0$.
En ese caso, nos queda el sistema $x \prime = - T_v$.
Tiene una solución trivial $x_{(t)} = x_{(0)} - T_v \cdot x$ y es estable no-asintóticamente.

Cuando $(T_c = 0, T_v \neq 0)$, la solución analítica es una recta con pendiente negativa.
En el momento $t_{crit} = x_{(0)} / \cdot T_v$, la solución cruza el eje de las abcisas.
Si el almacén es de un producto físico u otro almacenable que no pueda incurrir en deuda, y su stock no pueda ser nulo, de ahí en más $x = 0 \forall t \geq t_{crit}$.

Recordemos que el sistema original llega a converger asintóticamente cuando $T_c < T_v$.
Este aspecto de su comportamiento, como los otros que analizamos, son diferentes al comportamiento del sistema original.

\section{Conclusiones}

De los errores calculados podemos ordenar los métodos de la siguiente manera, de mayor a menor precisión:
\begin{itemize}
\item Runge-Kutta 4,5
\item Runge-Kutta 4
\item Runge-Kutta 2,3
\item Runge-Kutta 2
\item Euler
\end{itemize}

Las diferencias en los errores entre cada método son de entre uno y dos órdenes de magnitud y, como es de esperar, los métodos de mayor orden llevan errores menores, aproximando más cercanamente la solución analítica.

Podemos notar que Euler, siendo el mâs sencillo de todos, fue en tiempo significativamente más rápido que el resto.
Sin embargo, sus errores son varios órdenes mayores a los de los resultados  obtenidos con los métodos restantes, como el RK4.

Dada la baja complejidad computacional de los métodos implementados por nosotros al momento de correrlos, notamos que estos terminan antes su ejecución en comparación con los métodos estándar de mayor orden.
Sin embargo, los cálculos adicionales realizados en Runge-Kutta 2,3 y Runge-Kutta 4,5 como el uso de pasos variables demuestran sus resultados al darnos una mejor precisión en los resultados obtenidos.

Estas diferencias en la presición se deben al orden de cada método.
Podemos ver que los métodos de mayor orden aproximan mejor la solución. 
\end{document}
